class: center, middle, inverse, title-slide .title[ # Moving Beyond Linearity ] .subtitle[ ## Generalized Additive Models (GAM) ] .author[ ### Andreas Scharmüller ] .institute[ ### AG Landscape Ecology ] .date[ ### 2020-11-09 (updated: 2023-01-09) ] --- # Nonlinearity ```r ggplot(dat, aes(y = y, x = x)) + geom_point() ``` <!-- --> .foot-note[ <https://noamross.github.io/gams-in-r-course> ] --- class: center, inverse, middle # Linear Regression --- # Linear Model (LM) - Easy to interpret - Confined to linear relationships - Normally distributed responses `$$y_i = \beta_0 + \beta_1x_{1i} + \epsilon_i, \epsilon \sim N(0,\sigma^2)$$` --- # Linear Model (LM) .pull-left[ ```r lm(y ~ x, data = data) ``` ] .pull-right[ __R2:__ 0.17 __AIC:__ 513.9895637 ] <!-- --> --- # Generalized Linear Models (GLM) - Additional distributions (Poisson, Gamma, Binomial, etc.) `$$\mathbb{E}(y_i) = \beta_0 + \beta_1x_{1i} + \beta_2x_{2i} + \epsilon_i$$` --- # Generalized Linear Models (GLM) .pull-left[ ```r glm(y ~ x, data = data, family = 'Gamma') ``` ] .pull-right[ __R2:__ 0.13 __AIC:__ 531 ] <!-- --> --- # Polynomial Regression - Extra predictors, obtained by raising the original predictors to a power - Specific patterns, not very flexible - Global not Local - Might lead to poor residuals, predictions, extrapolations - Especially at the boundaries `$$y_i = \beta_0 + \beta_1x_{1i} + \beta_2x_{1i}^2 + \beta_3x_{1i}^3 + \epsilon_i$$` --- # Polynomial Regression .pull-left[ ```r lm(y ~ poly(x, 3), data = data) ``` ] .pull-right[ 3rd order | __R2:__ 0.5 __AIC:__ 467 9th order | __R2:__ 0.71 __AIC:__ 425 ] <!-- --> --- class: center, inverse, middle # Generalized Additive Models (GAMs) --- # Generalized Additive Models (GAMs) .pull-left[ ```r require(mgcv) gam(y ~ s(x), data = data) ``` ] .pull-right[ __Dev. explained:__ 0.69 __AIC:__ 425 ] <!-- --> --- # Generalized Additive Models (GAMs) `$$y_i = \beta_0 + f_1(x_{1i}) + f_2(x_{2i}) + f_3(x_{3i},x_{4i}) + ... + \epsilon_i, \epsilon \sim N(0,\sigma^2)$$` `$$\mathbb{E}(Y) = g^{-1} \bigg( \beta_0+ \sum_{j=1}^Jf_j(x_j)\bigg)$$` `$$f_j(x_j) = \sum_{k=1}^K \beta_{j,k} b_{j,k} (x_j)$$` - Smoothing function: `\(f()\)` - `\(b_{j,k}\)` - `\(K\)` Basis functions for each covariate `\(x_j\)` - `\(\beta_{j,k}\)` - `\(K\)` Coefficients for each covariate `\(x_j\)` --- # Generalized Additive Models (GAMs) <br><br><br> <img src="https://raw.githubusercontent.com/noamross/gams-in-r-course/master/images/tradeoff-slider.png"> .footnote[ <https://noamross.github.io/gams-in-r-course> ] --- # Generalized Additive Models (GAMs) <!-- NOTE doesn't work anymore --> <!-- tweetrmd::include_tweet("https://twitter.com/ucfagls/status/842444686513991680") --> <blockquote class="twitter-tweet" data-width="550" data-lang="en" data-dnt="true" data-theme="light"><p lang="en" dir="ltr">140 char vrsn<br><br>1 GAMs are just GLMs<br>2 GAMs fit wiggly terms<br>3 use + s(foo) not foo in frmla<br>4 use method = "REML"<br>5 gam.check()</p>&mdash; Dr Gavin Simpson 😷🇪🇺 (@ucfagls) <a href="https://twitter.com/ucfagls/status/842444686513991680?ref_src=twsrc%5Etfw">March 16, 2017</a></blockquote> .footnote[ <https://www.fromthebottomoftheheap.net> ] --- # Generalized Additive Models (GAMs) - Predictor replaced by `\(K\)` synthetic predictors (Basis Functions) - Basis Functions are Weighted by Coefficients - Summed to build Regression Line __(Spline)__ - Wiggliness is Penalized - Minimize Sums of Squared Error (SSE) - Minimize Wiggliness (Overfitting) `$$Fit = Likelihood - \lambda \times Wiggliness$$` `$$\hat{\beta} = arg\,min_\beta || y-X\beta||^2 + \lambda \beta^TS\beta$$` --- # Generalized Additive Models (GAMs) ### Example: Airquality ```r # data airqu = na.omit(airquality) head(airqu) ``` ``` ## Ozone Solar.R Wind Temp Month Day ## 1 41 190 7.4 67 5 1 ## 2 36 118 8.0 72 5 2 ## 3 12 149 12.6 74 5 3 ## 4 18 313 11.5 62 5 4 ## 7 23 299 8.6 65 5 7 ## 8 19 99 13.8 59 5 8 ``` ```r # model air1 = gam(Ozone ~ s(Wind, bs = 'cr', k = 7), # Cubic Regr. Spline data = airqu, method = 'REML') ``` --- # Generalized Additive Models (GAMs) ### Example: Airquality <!-- --> --- class: center, inverse, middle # Splines <img src="https://s3files.core77.com/blog/images/513221_81_55368_4VoVB83F4.jpg" width="50%"> .footnote[ <https://www.core77.com/posts/55368/When-Splines-Were-Physical-Objects> ] --- # Splines Predictor replaced by K synthetic predictors (Basis Functions) ```r predict(air1, type = 'lpmatrix') # or: model.matrix(air1) ``` <!-- --> --- # Splines Basis Functions are Weighted by Coefficients ```r predict(air1, type = 'lpmatrix') %*% diag(coef(air1)) ``` <!-- --> --- # Splines Summed to build Regression Line (Spline) <!-- --> --- # Splines ### Complexity ! ```r lm1 = lm(y ~ x, data = dat) coef(lm1) ``` ``` ## (Intercept) x ## 11.509559 -5.235207 ``` ```r gam1 = gam(y ~ s(x, k = 7), data = dat) coef(gam1) ``` ``` ## (Intercept) s(x).1 s(x).2 s(x).3 s(x).4 s(x).5 s(x).6 ## 8.885758 -8.004486 -12.522057 -5.455329 -3.854084 19.889658 6.737313 ``` --- # Splines ```r air2 = gam(Ozone ~ s(Wind, bs = 'cr', k = 3), # Cubic Regr. Spline data = airqu, method = 'REML') ``` ``` ## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0. ## ℹ Please use `linewidth` instead. ``` <!-- --> --- # Splines - Knots Check for k: ```r gam.check(air2) ``` ```r k.check(air1) # k = 7 ``` ``` ## k' edf k-index p-value ## s(Wind) 6 3.427874 0.8576004 0.0525 ``` ```r k.check(air2) # k = 3 ``` ``` ## k' edf k-index p-value ## s(Wind) 2 1.966923 0.8369715 0.0225 ``` -- If p-value is low, k might be too small --- # Overfitting? ### Large basis size (knots) lead to overfitting? ```r air3 = gam(Ozone ~ s(Wind, bs = 'cr', k = length(unique(airqu$Wind))), # 29 data = airqu, method = 'REML') ``` <!-- --> --- # Overfitting? SOLUTION: Penalize departure from smoothness: `$$P(f) = \int f''(x)^2 dx = \beta^T\mathbf{S}\beta$$` `$$Fit = Likelihood - \lambda \times Wiggliness$$` <!-- OLD formula --> <!-- `$$\sum_{i=1}^n (y_i - g(x_i))^2 + \lambda \int g''(t)^2dt$$` --> `$$\hat{\beta} = \operatorname{arg min}_\beta \|y-X\beta\|^2 + \lambda\beta^TS\beta$$` - Likelihood: How well a GAM captures patterns in the data - Wiggliness: Complexity of a smooth - Penalty matrix: `\(S\)` - Smoothing parameter `\(\lambda\)` is optimized in `gam()` - Controls the trade-off between Likelihood and Wiggliness .footnote[ <https://www.maths.ed.ac.uk/~swood34/mgcv/tampere/basis-penalty.pdf> ] --- # Smoothing ### Smoothing parameter ```r gam(y ~ s(x1, sp = NULL), data = data, sp = NULL, method = 'REML') # GCV.Cp, ML ``` - Smoothing parameter estimation method: `method =` - Recommended: `method = 'REML'`¹ .footnote[ [1] Wood et al. (2011) Maximum Likelihood Algorithm: <https://www.youtube.com/watch?v=XepXtl9YKwc> ] ??? GCV...Generalized Cross-Validation REML...Residual/Restricted Maximum Likelihood ML...Maximum Likelihood --- # Smoothing ### Smoothing parameter <!-- --> - A: `sp = 0` - B: `sp = 1e5` - C: `method = 'REML'` -- __The R-package does that for you!__ .footnote[ Smoothness selection alogrithms: <https://www.maths.ed.ac.uk/~swood34/mgcv/tampere/smoothness.pdf> ] --- class: center, inverse, middle # GAMs in R --- # GAMs in R ##### `gam` package - Original GAM R-package #### `mgcv` package - Most used package ##### `qgam` package - Quantile GAMs. Model e.g. the 75% percentile ##### `gamlss` package - Model not only the location & scale, but also the shape (e.g. kurtosis) ##### `brms` package - Bayesian GAM approach - runs [STAN](https://mc-stan.org/) --- # The mgcv package `mgcv::gam()` - Standard function `mgcv::bam()` - Reduces RAM-overhead `mgcv::gamm()` - Uses the `nlme` package for random effects `gamm4::gamm4()` - Uses the `lme4` package for random effects --- # GAMs in R R-function ```r require(mgcv) gam(y ~ s(x1, bs = 'tp', k = -1) + x2, # formula data = data, # data family = 'gaussian', # family object method = 'REML', # default: 'GCV.Cp' sp = NULL) # smoothing parameter ``` --- # GAMs in R Formula ```r require(mgcv) # gam(y ~ s(# smooth term s(), te() x1, # predictor bs = 'tp', # spline basis k = -1, # number of basis functions (i.e. knots) sp = NULL # smoothing parameter ) # + x2, # linear term # data = data, # family = 'gaussian', # method = 'REML' # default: 'GCV.Cp' # sp = NULL) ``` --- class: center, inverse, middle # Smooth terms --- # Smooth terms Two additive smooths ```r gam(y ~ s(x1) + s(x2)) ``` -- Smooth-interactions ```r gam(y ~ s(x1, x2)) # common way to declare spatial data gam(y ~ s(x1, by = fac)) ``` -- Tensor product smooths ```r gam(y ~ te(x1, x2, k = c(4,8))) # interaction on different scales ``` --- # Splines - Thin Plate Regression Splines (TPRS) - Default - Computationally more demanding (CRS) - Cubic Regression Splines (CRS) - Computationally less demanding - Cyclic Cubic Regression Splines - Random Effects - Different penalty matrices: `\(S\)` .footnote[ <https://stats.stackexchange.com/questions/305338/adaptive-gam-smooths-in-mgcv> ] <!-- --- --> <!-- # Splines --> <!-- __Thin Plate__ Regression Splines --> <!-- - Any number of dimensions --> <!-- - Smoothing is the same (isotropy, rotational invariance) --> <!-- - Second-order: --> <!-- - Penalty is porportionate to the integral of the squared second derivative --> <!-- - `\(m=1\)`: Penalizing squared first derivatives. --> ```r gam(y ~ s(x, bs = 'tp'), data = data) ``` --- # Splines __Cyclic Cubic__ regression splines - Cyclical data (e.g. seasons) ```r gam(y ~ s(x, bs = 'cc'), data = data) ``` <img src="data/pre_emergence_herbicides.png" width="500px" style="position:absolute; right:30px; bottom:20px;"> --- # Splines __Soap Films__ - Boundary polygons can be introduced - Spatial models ```r gam(y ~ s(x, y, bs = 'so', xt = list(bnd = boundary_polygon)), data = data) ``` <figure> <img src="https://fromthebottomoftheheap.net/assets/img/posts/soap-film-smoothers-ggplot-soap-film-1.png" width="80%"> <caption>https://fromthebottomoftheheap.net/2016/03/27/soap-film-smoothers</caption> </figure> --- class: center, inverse, middle # Random Effects --- # Splines Discrete __random effects__ - Classes (e.g. age, sex) - Sites, states, rivers, lakes - Mesocosm groups - No need to set k (equals number of levels) ```r gam(y ~ s(x) + s(fac, bs = 're'), data = data) ``` --- # gamm() & gamm4() - Calls `nlme::lme()` ```r mgcv::gamm(y ~ s(x), data = data, random = list(fac = ~1)) ``` - Calls `lme4::lmer()` ```r gamm4::gamm4(y ~ s(x), data = data, random = list(fac = ~ 1 | fac)) ``` --- # Random Effects "Random effects and penalized splines are the same thing." <!-- tweetrmd::include_tweet("https://twitter.com/ericJpedersen/status/1293508069016637440") --> <!-- --> .footnote[ <https://twitter.com/ericJpedersen/status/1293508069016637440> <https://www.tjmahr.com/random-effects-penalized-splines-same-thing> ] --- class: center, inverse, middle # Model example --- # Model example ```r require(gamair) data('mack') head(mack, n = 3) ``` ``` ## egg.count egg.dens b.depth lat lon time salinity flow s.depth temp.surf temp.20m net.area ## 1 0 0 4342 44.57 -4.65 8.23 35.72 417 104 15.0 15.0 0.242 ## 2 0 0 4334 44.57 -4.48 9.68 35.74 405 98 15.4 15.4 0.242 ## 3 0 0 4286 44.57 -4.30 10.90 35.74 377 101 15.9 15.9 0.242 ## country vessel vessel.haul c.dist ## 1 SP COSA 22 0.8395141 ## 2 SP COSA 23 0.8591926 ## 3 SP COSA 24 0.8930153 ``` -- - Response: `egg.count` - Covariates: - Sea bed depth at sampling location: `b.depth` - Water salinity: `salinity` .footnote[ <https://cran.r-project.org/web/packages/gamair/gamair.pdf> ] --- # Model example ```r mod = gam(egg.count ~ s(b.depth) + s(salinity), data = mack, family = 'nb') ``` --- # Model example Summary ```r summary(mod) ``` Predict ```r predict(mod) ``` AIC ```r AIC(mod) ``` Checking ```r gam.check(mod) k.check(mod) ``` --- # Model summary ``` ## ## Family: Negative Binomial(0.455) ## Link function: log ## ## Formula: ## egg.count ~ s(b.depth) + s(salinity) ## ## Parametric coefficients: ## Estimate Std. Error z value Pr(>|z|) ## (Intercept) 2.03383 0.08967 22.68 <2e-16 *** ## --- ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 ## ## Approximate significance of smooth terms: ## edf Ref.df Chi.sq p-value ## s(b.depth) 7.045 8.076 99.21 <2e-16 *** ## s(salinity) 6.672 7.740 90.90 <2e-16 *** ## --- ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 ## ## R-sq.(adj) = 0.151 Deviance explained = 37.4% ## -REML = 1047.9 Scale est. = 1 n = 330 ``` --- # Model prediction ```r # newdata new = with(mack, expand.grid(b.depth = mean(b.depth), salinity = seq(min(salinity, na.rm = TRUE), max(salinity, na.rm = TRUE), length.out = 100))) require(data.table) setDT(new) # convert to data.table # predict prd = predict(mod, newdata = new, type = 'response', exclude = 's(b.depth)', # exclude of a variable se.fit = TRUE) # include standard errors # update newdata new[ , `:=` (fit = prd$fit, lwr = prd$fit - (1.96 * prd$se.fit), upr = prd$fit + (1.96 * prd$se.fit)) ] ``` --- # Model plot ```r ggplot(new, aes(y = fit, x = salinity)) + geom_point(data = mack, aes(y = egg.count)) + geom_ribbon(aes(ymin = lwr, ymax = upr), alpha = 0.4) + geom_line(size = 1.5, col = 'blue') ``` <!-- --> --- # Model checking <!-- --> ``` ## ## Method: REML Optimizer: outer newton ## full convergence after 6 iterations. ## Gradient range [-1.05935e-06,1.198944e-06] ## (score 1047.905 & scale 1). ## Hessian positive definite, eigenvalue range [0.6362334,106.8741]. ## Model rank = 19 / 19 ## ## Basis dimension (k) checking results. Low p-value (k-index<1) may ## indicate that k is too low, especially if edf is close to k'. ## ## k' edf k-index p-value ## s(b.depth) 9.00 7.04 0.78 0.18 ## s(salinity) 9.00 6.67 0.60 <2e-16 *** ## --- ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 ``` ??? gam.check() - console output - number of basis functions - If a term in gam.check is significant, increas k, refit and check again gam.check() - plots - Q-Q plot (compares residuals to a normal distribution) - follow straight line - Histogram of residuals - should be bell symmetrically shaped - Residual vs. Linear predictors - Should be evenly distributed around zero - Response vs. Fitted values - Perfect model would be a straight line - We would also be happy if it would follow a 1-to-1 line --- # Model checking <!-- --> ``` ## ## Method: REML Optimizer: outer newton ## full convergence after 5 iterations. ## Gradient range [-0.0001075229,9.593977e-05] ## (score 71.222 & scale 4.071951). ## Hessian positive definite, eigenvalue range [0.04951955,14.74743]. ## Model rank = 19 / 19 ## ## Basis dimension (k) checking results. Low p-value (k-index<1) may ## indicate that k is too low, especially if edf is close to k'. ## ## k' edf k-index p-value ## s(hp) 9.00 1.36 0.95 0.37 ## s(disp) 9.00 4.67 1.37 0.95 ``` --- # Model selection - Adds additional shrinkage on perfectly smooth terms - Shrinkage Smoothers - `bs = 'ts'`, `bs = 'cs'` - Double penalty approach - `select = TRUE` - Adds second penalty to all terms - Penalizes the NULL space - "Shrinkage" (Ridge Regression, Lasso) ```r gam(y ~ s(x1) + s(x2, bs = 'ts') + s(x3), select = TRUE, data = data, method = 'REML') ``` .footnote[ Marra & Wood (2011) ] --- # Model selection - AIC - Expert judgement - Computational time - __Inferential goals of the study__ .footnote[ <https://stats.stackexchange.com/questions/274151/anova-to-compare-models> ] --- class: center, middle # Exercise 1 ## GAM/exercise_pesticides.Rmd --- class: center, inverse, middle # Visualization --- # Visualization ```r ir1 = gam(Sepal.Length ~ s(Petal.Length), data = iris, family = 'gaussian', method = 'REML') ``` --- # Visualization Partial effect plots .pull-left[ ```r require(mgcv) plot(ir1, pages = 1) ``` <!-- --> ] .pull-right[ ```r require(gratia) draw(ir1) ``` <!-- --> ] --- # Visualization Multiple covariates ```r mt1 = gam(mpg ~ s(disp) + s(hp), data = mtcars, family = 'gaussian', method = 'REML') ``` --- # Visualization ```r vis.gam(mt1, # GAM object view = c("disp", "hp"), # variables plot.type = "persp", # 3D plot theta = 135, # horizontal rotation phi = 10, # phi vertical rotation r = 10) # zoom ``` <!-- --> ??? yellow: larger predictions<br> red: smaller predictions<br> mpg: miles per gallon<br> disp: displacement (German: Hubraum)<br> hp: horsepower<br> --- # Visualization ```r vis.gam(mt1, # GAM object view = c("disp", "hp"), # variables plot.type = "contour") # contour plot or heatmap ``` <!-- --> --- class: center, inverse, middle # Hierachical / Mixed-effect # GAMs --- # Hierarchical GAMs - Hierarchical data: grouped data (factor!) - Non-linear relationships - Global function - Group-specific function <br>  --- # Hierarchical GAMs <figure> <img src="data/hgam_concept.png" width="80%"> <caption>Pedersen et al. (2019)</caption> </figure> - Estimating a function for each species throws away shared information - Highly noisy estimates - Ignoring inherent grouping would miss individual optima --- # Hierarchical GAMs - Should there be a common smoother or a smoother for each group? - Do all group-specific smoothers have the __same wiggliness__ or should each group wiggly independently (own smoothing parameter)? - Will the smoothers for each group have a __similar shape__? --- # Hierarchical GAMs <figure> <img src="data/hgam_possibilities.png" width="60%"> <caption>Pedersen et al. (2019)</caption> </figure> --- ## Hierarchical GAM (HGAM) ### Chlortoluron ```r dat[ , head(.SD, 2), site_id] ``` ``` ## site_id value_fin date name doy site_id_f ## <char> <num> <Date> <char> <int> <fctr> ## 1: RP_2619521200 0.110 2011-11-28 Chlortoluron 332 RP_2619521200 ## 2: RP_2619521200 0.310 2011-12-12 Chlortoluron 346 RP_2619521200 ## 3: RP_2649525000 0.052 2010-02-22 Chlortoluron 53 RP_2649525000 ## 4: RP_2649525000 0.140 2010-02-08 Chlortoluron 39 RP_2649525000 ## 5: SN_OBF00200 0.068 2001-01-31 Chlortoluron 31 SN_OBF00200 ## 6: SN_OBF00200 0.053 2001-03-01 Chlortoluron 60 SN_OBF00200 ## 7: SN_OBA02810 0.022 2010-04-25 Chlortoluron 115 SN_OBA02810 ## 8: SN_OBA02810 0.012 2010-04-07 Chlortoluron 97 SN_OBA02810 ``` --- ## Hierarchical GAM (HGAM) ### Chlortoluron <!-- --> --- ## No HGAM ```r m0 = gam(value_fin ~ s(doy, bs = 'cc', k = 12), data = dat, family = Gamma(link = 'log')) ``` <!-- --> --- ## HGAM - Global Smooth ```r mG = gam(value_fin ~ s(doy, bs = 'cc', k = 12) + s(site_id_f, bs = 're'), data = dat, family = Gamma(link = 'log')) ``` - Single global smooth term for each variable - Level-individual __random effect intercepts__ (`bs = 're'`) - Mixed models (like in lme4, nlme) --- ## HGAM - Global Smooth ```r mG = gam(value_fin ~ s(doy, bs = 'cc', k = 12) + s(site_id_f, bs = 're'), data = dat, family = Gamma(link = 'log')) ``` <!-- --> --- ## HGAM - Global + Shared Group Level Smooth ```r mGS = gam(value_fin ~ s(doy, bs = 'cc', k = 12) + s(doy, site_id_f, bs = 'fs', k = 12, xt=list(bs='cc')) + s(site_id_f, bs = 're'), # not explicitly needed data = dat, family = Gamma(link = 'log')) ``` - Factor-smooth interaction (`bs = 'fs'`) - Single global smooth term for each variable - Varying slopes - Copy of each set of basis functions, - Penalizing functions too far away from the average - __Similar wiggliness__, one smoothing parameter --- ## HGAM - Global + Shared Group Level Smooth ```r mGS = gam(value_fin ~ s(doy, bs = 'cc', k = 12) + s(doy, site_id_f, bs = 'fs', k = 12, xt=list(bs='cc')) + s(site_id_f, bs = 're'), # not explicitly needed data = dat, family = Gamma(link = 'log')) ``` <!-- --> --- ## HGAM - Global + Individual Group Level Smooth ```r mGI = gam(value_fin ~ s(doy, bs = 'cc', k = 12) + s(doy, by = site_id_f, bs = 'cc', k = 12) + s(site_id_f, bs = 're'), # explicitly needed data = dat, family = Gamma(link = 'log')) ``` - `s(x, by = fac)` and `s(fac, bs = 're')` - Similar to GS - Varying slopes - __Different wiggliness__, multiple smoothing parameters - Useful: Different wiggliness expected --- ## HGAM - Global + Individual Group Level Smooth ```r mGI = gam(value_fin ~ s(doy, bs = 'cc', k = 12) + s(doy, by = site_id_f, bs = 'cc', k = 12) + s(site_id_f, bs = 're'), # explicitly needed data = dat, family = Gamma(link = 'log')) ``` <!-- --> --- ## HGAM - Only Shared Group Level Smooth ```r mS = gam(value_fin ~ s(doy, site_id_f, bs = 'fs', k = 12, xt=list(bs='cc')) + s(site_id_f, bs = 're'), # not explicitly needed data = dat, family = Gamma(link = 'log')) ``` <!-- --> --- ## HGAM - Only Individual Group Level Smooth ```r mI = gam(value_fin ~ s(doy, by = site_id_f, bs = 'cc', k = 12) + s(site_id_f, bs = 're'), # explicitly needed data = dat, family = Gamma(link = 'log')) ``` <!-- --> --- class: center, middle # Exercise 2 ## GAM/exercise_co2uptake.Rmd --- # Tutorials, Blogs & Videos - Generalized Additive Models in R (Noam Ross) - <https://noamross.github.io/gams-in-r-course> - Doing magic and analyzing seasonal time series with GAM (Generalized Additive Model) in R - <https://petolau.github.io/Analyzing-double-seasonal-time-series-with-GAM-in-R/> - From the Bottom of the Heap - Blog by Gavin Simpson - <https://www.fromthebottomoftheheap.net/> - Introduction lecture - By Gavin Simpson - <https://www.youtube.com/watch?v=sgw4cu8hrZM> - Introducing gratia - <https://www.fromthebottomoftheheap.net/2018/10/23/introducing-gratia/> - Noam Ross Github - <https://github.com/noamross/gam-resources> - Environmental computing - <http://environmentalcomputing.net/intro-to-gams> - Random effects and penalized splines are the same thing - <https://www.tjmahr.com/random-effects-penalized-splines-same-thing> - GAM course - <https://noamross.github.io/mgcv-esa-2018> --- # Slides - OLAT - <https://andschar.github.io/teaching> ## Made with - <https://github.com/rstudio/rmarkdown> - <https://github.com/yihui/knitr> - <https://github.com/yihui/xaringan> --- # References Pedersen et al. (2019) Marra & Wood (2011) Wood (2011) Burnham & Anderson (1998) --- layout: false # [Moving Beyond Linearity](https://andschar.github.io/teaching/GAM-intro.html) ## Thank you for your attention! <https://andschar.github.io/teaching> __Andreas Scharmüller__<br> Quantitative Landscape Ecology<br/> iES Landau, Institute for Environmental Sciences<br> University of Koblenz-Landau<br> University of Strasbourg <img src="data:image/png;base64,#https://upload.wikimedia.org/wikipedia/commons/thumb/4/4f/Twitter-logo.svg/500px-Twitter-logo.svg.png" width="20"> @andschar <br/> <img src="https://seeklogo.com/images/M/mail-icon-logo-28FE0635D0-seeklogo.com.png" width="20"> andschar@proton.me <img src="https://www.uni-koblenz-landau.de/de/landau/fb7/umweltwissenschaften/logos/land-ecol-neu.jpg" width="100"> <img src="https://www.uni-koblenz-landau.de/de/uni/organisation/verwaltung/abteilungen/abt-1/dokumente-ab1/uni-logo-farbig.jpg" width="200"> <!-- OLD MATERIAL --> <!-- --- --> <!-- # Mixed Models --> <!-- $$ --> <!-- Y = X \beta + Zb + \epsilon \\ --> <!-- b \sim Normal(0, \sigma_b) \\ --> <!-- \epsilon \sim Normal(0, \sigma_y) \\ --> <!-- X: \textrm{fixed effects model matrix} \\ --> <!-- Y: \textrm{mixed effects model matrix} \\ --> <!-- \sigma_b, \sigma_y: \textrm{variance components} \\ --> <!-- \sigma_b: \textrm{where the magic happens} --> <!-- $$ --> <!-- From: <https://www.tjmahr.com/random-effects-penalized-splines-same-thing> --> <!-- --- --> <!-- # Basis functions --> <!-- ```{r eval=FALSE} --> <!-- gam(y ~ s(x1, k = 3), --> <!-- data = data) --> <!-- ``` --> <!-- ```{r echo=FALSE, fig.width=6, fig.height=3} --> <!-- oc1 = gam(SP ~ s(P, k = 3), --> <!-- data = oc, --> <!-- method = 'REML') --> <!-- oc2 = gam(SP ~ s(P, k = -1), --> <!-- data = oc, --> <!-- method = 'REML') --> <!-- oc$pr1 = predict(oc1, type = 'response') --> <!-- oc$pr2 = predict(oc2, type = 'response') --> <!-- ggplot(oc, aes(x = P)) + --> <!-- geom_point(aes(y = SP)) + --> <!-- geom_line(aes(y = pr1), col = 'blue') + --> <!-- labs(x = 'x', y = 'y') --> <!-- ``` --> <!-- --- --> <!-- # Basis functions --> <!-- ```{r eval=FALSE} --> <!-- gam(y ~ s(x1), --> <!-- data = data) --> <!-- ``` --> <!-- ```{r echo=FALSE, fig.width=6, fig.height=3} --> <!-- oc1 = gam(SP ~ s(P, k = 3), --> <!-- data = oc, --> <!-- method = 'REML') --> <!-- oc2 = gam(SP ~ s(P, k = -1), --> <!-- data = oc, --> <!-- method = 'REML') --> <!-- oc$pr1 = predict(oc1, type = 'response') --> <!-- oc$pr2 = predict(oc2, type = 'response') --> <!-- ggplot(oc, aes(x = P)) + --> <!-- geom_point(aes(y = SP)) + --> <!-- geom_line(aes(y = pr1), col = 'blue') + --> <!-- geom_line(aes(y = pr2), col = 'red') + --> <!-- labs(x = 'x', y = 'y') --> <!-- ``` --> <!-- --- --> <!-- # Global Smoother --> <!-- - Single global smooth term for each variable --> <!-- - Level-individual __random effect intercepts__ (`bs = 're'`) --> <!-- - Mixed models (like in lme4, nlme) --> <!-- ```{r} --> <!-- spr_G = gam(logSize ~ s(days) + --> <!-- s(plot, bs = 're'), --> <!-- data = Spruce, --> <!-- method = 'REML') --> <!-- Spruce$fit_G = fitted(spr_G) --> <!-- ggplot(Spruce) + --> <!-- geom_line(aes(y = fit_G, x = days, col = plot)) + --> <!-- labs(y = 'Fit') --> <!-- ``` --> <!-- ```{r, echo=FALSE, warning=FALSE, fig.width=5, fig.height=2.5} --> <!-- plot(spr_G, pages = 1, main = NULL) --> <!-- ``` --> <!-- --- --> <!-- # Global Smoother + Group Level Smoother --> <!-- - Factor-smooth interaction (`bs = 'fs'`) --> <!-- - Single global smooth term for each variable --> <!-- - __Varying slopes__ --> <!-- - Same wiggliness (i.e. complexity) of smooths --> <!-- ```{r} --> <!-- spr_GS = gam(logSize ~ s(days, plot, bs = 'fs'), --> <!-- data = Spruce, --> <!-- method = 'REML') --> <!-- ``` --> <!-- ```{r, echo=FALSE, warning=FALSE, fig.width=5, fig.height=2.5} --> <!-- plot(spr_GS, pages = 1, main = NULL) --> <!-- ``` --> <!-- --- --> <!-- # Global Smoother + Group Level Smoother --> <!-- - `s(x, by = fac)` and `s(fac, bs = 're')` --> <!-- - Similar to GS --> <!-- - Varying slopes --> <!-- - __Different wiggliness__ (i.e. complexity) of smooths --> <!-- ```{r} --> <!-- spr_GI = gam(logSize ~ s(days) + --> <!-- s(days, by = plot) + --> <!-- s(plot, bs = 're'), --> <!-- data = Spruce, --> <!-- method = 'REML') --> <!-- ``` --> <!-- ```{r, echo=FALSE, warning=FALSE, fig.width=5, fig.height=3.5} --> <!-- plot(spr_GI, pages = 1) --> <!-- ``` --> <!-- --- --> <!-- # Splines --> <!-- ```{r echo=FALSE, fig.height=5} --> <!-- # TODO how should the knots be ordered? --> <!-- k = 7 --> <!-- knots = data.frame(x = seq(min(dat$x), max(dat$x), length.out = k)) --> <!-- spl = smoothCon(s(x, bs = 'cr', k = k), --> <!-- data = dat,#[ ,'x', drop = FALSE], --> <!-- knots = knots)[[1]] --> <!-- coefs = coef(gam(y ~ s(x, bs = 'cr', k = k), data = dat)) --> <!-- fx = spl$X %*% coefs --> <!-- bs_coef = sweep(spl$X, 2, coefs, "*") --> <!-- bs1 = sweep(spl$X, 2, rep(1, k), "*") --> <!-- # plot --> <!-- ylim <- range(fx, bs_coef) --> <!-- op <- par(mar = c(5,4,1,2) + 0.1) --> <!-- layout(matrix(1:2, ncol = 2)) --> <!-- # 1 --> <!-- plot(fx ~ dat$x, type = "n", ylim = ylim, ylab = "", xlab = "x") --> <!-- abline(v = knots$x, lty = "dotted") --> <!-- matlines(dat$x, bs1, lwd = 2) --> <!-- title(ylab = expression(f(x)), line = 2.75) --> <!-- # 2 --> <!-- plot(fx ~ dat$x, type = "n", ylim = ylim, ylab = "", xlab = "x") --> <!-- abline(v = knots$x, lty = "dotted") --> <!-- matlines(dat$x, bs_coef, lwd = 2) --> <!-- lines(fx ~ dat$x, lwd = 4) --> <!-- title(ylab = expression(f(x)), line = 2.75) --> <!-- par(op) --> <!-- layout(1) --> <!-- ``` --> <!-- --- --> <!-- # Splines --> <!-- - Splines replace a predictor variable, with a set of synthetic predictor variables. --> <!-- - Family of functions or transformations that can be applied to `\(x\)` --> <!-- - Polynomial basis function --> <!-- `$$f_j(x_i) = x_i^j$$` --> <!-- - Piecewise constant basis function: --> <!-- `$$f_j(x_i) = I(c_j \leq x_i \leq c_{j+1})$$` --> <!-- - Penalized splines: reduce wiggliness --> <!-- .footnote[ --> <!-- https://www.youtube.com/watch?v=ENxTrFf9a7c&t=2226 --> <!-- https://www.tjmahr.com/random-effects-penalized-splines-same-thing/ --> <!-- ] -->